{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Petroleum Blending (Solved)\n", "\n", "## Try me\n", " [![Open In Colab](../../_static/colabs_badge.png)](https://colab.research.google.com/github/ffraile/operations-research-notebooks/blob/main/docs/source/CLP/solved/Petroleum%20Blending%20(Solved%20CBC).ipynb)[![Binder](../../_static/binder_badge.png)](https://mybinder.org/v2/gh/ffraile/operations-research-notebooks/main?labpath=docs%2Fsource%2FCLP%2Fsolved%2FPetroleum%20Blending%20(Solved%20CBC).ipynb)\n", "\n", "## Problem Definition\n", "Harry Stamper is an engineer working on a petroleum company. The company produces three grades of motor oil – super, premium, and extra. The three grades are made from the same 3 ingredients, named component 1, component 2 and component 3. The company wants to determine the optimal mix of the 3 components in each grade of motor oil that will maximize profit. \n", "\n", "The maximum availability of each component and their cost per barrel are as follows:\n", "\n", "| Component | Max barrels available / day | Cost / barrel |\n", "|--------------|------------------------------|---------------|\n", "| Component 1 | 4.500 | 12€ |\n", "| Component 2 | 2.700 | 10€ |\n", "| Component 3 | 3.500 | 14€ |\n", "\n", "To ensure the appropriate blend, each grade must have a minimum amount of component 1 plus a combination of other components as follows:\n", "\n", "| Grade | Component 1 | Component 2 | Component 3 | Selling price/barrel |\n", "|---------|--------------|------------------|-------------------|-----------------------|\n", "| Super | At least 50% | No more than 30% | - | 23€ |\n", "| Premium | At least 40% | - | No more than 25% | 20€ |\n", "| Extra | At least 60% | At least 10% | - | 18€ |\n", "\n", "The company wants to make at least 3000 barrels of each blend.\n", "*Formulate a LP to help Harry Stamper find the optimal blend that maximises profit.*" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Problem Model\n", "The decision variables are: \n", "\n", "$x_{ij}$ = amount in barrels of component i in grade j\n", "\n", "Objective function\n", "\n", "max $z = 23*(x_{11} + x_{21} + x_{31}) + 20*(x_{12} + x_{22} + x_{32})+18*(x_{13} + x_{23} + x_{33})-12*(x_{11} + x_{12} + x_{13}) - 10*(x_{21} + x_{22} + x_{23})-14*(x_{31} + x_{32} + x_{33})$\n", "\n", "Constraints\n", "\n", "$x_{11} \\geq 0.5*(x_{11} + x_{21} + x_{31})$ Component 1 requirement in Super\n", "\n", "$x_{21} \\leq 0.3*(x_{11} + x_{21} + x_{31})$ Component 2 requirement in Super\n", "\n", "$x_{12} \\geq 0.4*(x_{12} + x_{22} + x_{32})$ Component 1 requirement in Premium\n", "\n", "$x_{32} \\leq 0.25*(x_{12} + x_{22} + x_{32})$ Component 3 requirement in Premium\n", "\n", "$x_{13} \\geq 0.6*(x_{13} + x_{23} + x_{33})$ Component 1 requirement in Extra\n", "\n", "$x_{23} \\geq 0.1*(x_{13} + x_{23} + x_{33})$ Component 2 requirement in Extra\n", "\n", "\n", "$x_{11}+x_{12}+x_{13} \\leq 4500$ Component 1 availability\n", "\n", "$x_{21}+x_{22}+x_{23} \\leq 2700$ Component 1 availability\n", "\n", "$x_{31}+x_{32}+x_{33} \\leq 3500$ Component 1 availability\n", "\n", "$x_{11}+x_{21}+x_{31} \\geq 3000$ Super minimum production\n", "\n", "$x_{12}+x_{22}+x_{32} \\geq 3000$ Premium minimum production\n", "\n", "$x_{13}+x_{23}+x_{33} \\geq 3000$ Extra minimum production" ] }, { "cell_type": "code", "execution_count": null, "outputs": [], "source": [ "!pip install pandas\n", "!pip install numpy\n", "!pip install pulp" ], "metadata": { "collapsed": false, "pycharm": { "name": "#%%\n" } } }, { "cell_type": "code", "execution_count": 50, "metadata": {}, "outputs": [], "source": [ "import pandas as pd\n", "import numpy as np\n", "from IPython.display import display, Markdown\n", "import pulp" ] }, { "cell_type": "code", "execution_count": 51, "metadata": {}, "outputs": [], "source": [ "# Instantiate our problem class\n", "model = pulp.LpProblem(\"Maximising profits for Harry Stamper\", pulp.LpMaximize)" ] }, { "cell_type": "code", "execution_count": 52, "metadata": {}, "outputs": [], "source": [ "# Construct our decision variable lists\n", "components = ['Component 1', 'Component 2', 'Component 3']\n", "blends = ['Super', 'Premium', 'Extra']\n", "#decision variables\n", "variables = pulp.LpVariable.dicts(\"Quantity in barrels\",\n", " ((i, j) for i in blends for j in components),\n", " lowBound=0,\n", " cat='Continuous')\n", "\n", "# barrels is Super Component 1, super component 2, ..., extra component 3: x11, x21, x31, ...\n", "cost = [12, 10, 14]\n", "price = [23, 20, 18]\n", "coefficients =[price[i] - cost[j] for i in range(len(price)) for j in range(len(cost))]\n", "# Objective Function\n", "model += (\n", " pulp.lpSum([\n", " ((price[i]-cost[j]) * variables[(blends[i], components[j])]\n", " for i in range(len(blends)) for j in range(len(components)))])\n", ")\n", "\n", "# Requirements\n", "\n", "model += variables['Super','Component 1']>= 0.5*pulp.lpSum([variables['Super',components[j]] for j in range(len(components))]), \"Component 1 requirement in super\"\n", "model += variables['Super','Component 2']<= 0.3*pulp.lpSum([variables['Super',components[j]] for j in range(len(components))]), \"Component 2 requirement in super\"\n", "\n", "model += variables['Premium','Component 1']>= 0.4*pulp.lpSum([variables['Premium',components[j]] for j in range(len(components))]), \"Component 1 requirement in premium\"\n", "model += variables['Premium','Component 3']<= 0.25*pulp.lpSum([variables['Premium',components[j]] for j in range(len(components))]), \"Component 3 requirement in premium\"\n", "\n", "model += variables['Extra','Component 1']>= 0.6*pulp.lpSum([variables['Extra',components[j]] for j in range(len(components))]), \"Component 1 requirement in Extra\"\n", "model += variables['Extra','Component 2']>= 0.1*pulp.lpSum([variables['Extra',components[j]] for j in range(len(components))]), \"Component 2 requirement in Extra\"\n", "\n", "\n", "model += pulp.lpSum([variables[blends[i], 'Component 1'] for i in range(len(blends))]) <= 4500, \"Component 1 availability\" \n", "model += pulp.lpSum([variables[blends[i], 'Component 2'] for i in range(len(blends))]) <= 2700, \"Component 2 availability\" \n", "model += pulp.lpSum([variables[blends[i], 'Component 3'] for i in range(len(blends))]) <= 3500, \"Component 3 availability\"\n", "\n", "model += pulp.lpSum([variables['Super', components[j]] for j in range(len(components))]) >= 3000, \"Super minimum production\" \n", "model += pulp.lpSum([variables['Premium', components[j]] for j in range(len(components))]) >= 3000, \"Premium minimum production\" \n", "model += pulp.lpSum([variables['Extra', components[j]] for j in range(len(components))]) >= 3000, \"Extra minimum production\" \n" ] }, { "cell_type": "code", "execution_count": 53, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "'Optimal'" ] }, "execution_count": 53, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Solve our problem\n", "model.solve(solver=pulp.solvers.GUROBI(msg = 0))\n", "pulp.LpStatus[model.status]" ] }, { "cell_type": "code", "execution_count": 54, "metadata": { "pycharm": { "name": "#%%\n" } }, "outputs": [ { "data": { "text/markdown": [ "Total profit is 76800.00 €" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/markdown": [ "The following table shows the decision variables: " ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
VariablesSolution (GRB)Reduced cost (GRB)Objective Coefficient (GRB)Objective Lower bound (GRB)Objective Upper bound (GRB)
(Super, Component 1)Quantity_in_barrels_('Super',_'Component_1')1500.0000.00011.0008.000Inf
(Super, Component 2)Quantity_in_barrels_('Super',_'Component_2')450.0000.00013.00013.00013.000
(Super, Component 3)Quantity_in_barrels_('Super',_'Component_3')1050.0000.0009.0009.0009.000
(Premium, Component 1)Quantity_in_barrels_('Premium',_'Component_1')1200.0000.0008.000-Inf11.000
(Premium, Component 2)Quantity_in_barrels_('Premium',_'Component_2')1050.0000.00010.000-Inf10.000
(Premium, Component 3)Quantity_in_barrels_('Premium',_'Component_3')750.0000.0006.0006.00010.800
(Extra, Component 1)Quantity_in_barrels_('Extra',_'Component_1')1800.0000.0006.000-Inf17.333
(Extra, Component 2)Quantity_in_barrels_('Extra',_'Component_2')1200.0000.0008.0008.00025.000
(Extra, Component 3)Quantity_in_barrels_('Extra',_'Component_3')0.000-0.0004.000-Inf4.000
\n", "
" ], "text/plain": [ " Variables \\\n", "(Super, Component 1) Quantity_in_barrels_('Super',_'Component_1') \n", "(Super, Component 2) Quantity_in_barrels_('Super',_'Component_2') \n", "(Super, Component 3) Quantity_in_barrels_('Super',_'Component_3') \n", "(Premium, Component 1) Quantity_in_barrels_('Premium',_'Component_1') \n", "(Premium, Component 2) Quantity_in_barrels_('Premium',_'Component_2') \n", "(Premium, Component 3) Quantity_in_barrels_('Premium',_'Component_3') \n", "(Extra, Component 1) Quantity_in_barrels_('Extra',_'Component_1') \n", "(Extra, Component 2) Quantity_in_barrels_('Extra',_'Component_2') \n", "(Extra, Component 3) Quantity_in_barrels_('Extra',_'Component_3') \n", "\n", " Solution (GRB) Reduced cost (GRB) \\\n", "(Super, Component 1) 1500.000 0.000 \n", "(Super, Component 2) 450.000 0.000 \n", "(Super, Component 3) 1050.000 0.000 \n", "(Premium, Component 1) 1200.000 0.000 \n", "(Premium, Component 2) 1050.000 0.000 \n", "(Premium, Component 3) 750.000 0.000 \n", "(Extra, Component 1) 1800.000 0.000 \n", "(Extra, Component 2) 1200.000 0.000 \n", "(Extra, Component 3) 0.000 -0.000 \n", "\n", " Objective Coefficient (GRB) \\\n", "(Super, Component 1) 11.000 \n", "(Super, Component 2) 13.000 \n", "(Super, Component 3) 9.000 \n", "(Premium, Component 1) 8.000 \n", "(Premium, Component 2) 10.000 \n", "(Premium, Component 3) 6.000 \n", "(Extra, Component 1) 6.000 \n", "(Extra, Component 2) 8.000 \n", "(Extra, Component 3) 4.000 \n", "\n", " Objective Lower bound (GRB) Objective Upper bound (GRB) \n", "(Super, Component 1) 8.000 Inf \n", "(Super, Component 2) 13.000 13.000 \n", "(Super, Component 3) 9.000 9.000 \n", "(Premium, Component 1) -Inf 11.000 \n", "(Premium, Component 2) -Inf 10.000 \n", "(Premium, Component 3) 6.000 10.800 \n", "(Extra, Component 1) -Inf 17.333 \n", "(Extra, Component 2) 8.000 25.000 \n", "(Extra, Component 3) -Inf 4.000 " ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/markdown": [ "The following table shows the constraints: " ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
ConstraintRight Hand SideShadow PriceSlackMin RHSMax RHS
0Component_1_requirement_in_super0.00-18.000.00-850.000.00
1Component_2_requirement_in_super0.000.00450.00-450.00Inf
2Component_1_requirement_in_premium0.00-18.000.00-450.000.00
3Component_3_requirement_in_premium0.000.000.00-450.00450.00
4Component_1_requirement_in_Extra0.00-18.000.00-450.000.00
5Component_2_requirement_in_Extra0.000.00-900.00-Inf900.00
6Component_1_availability4500.0020.000.004500.006200.00
7Component_2_availability2700.004.000.002250.003150.00
8Component_3_availability3500.000.001700.001800.00Inf
9Super_minimum_production3000.000.00-0.00-Inf3000.00
10Premium_minimum_production3000.00-1.200.00-0.003000.00
11Extra_minimum_production3000.00-6.800.000.003000.00
\n", "
" ], "text/plain": [ " Constraint Right Hand Side Shadow Price Slack \\\n", "0 Component_1_requirement_in_super 0.00 -18.00 0.00 \n", "1 Component_2_requirement_in_super 0.00 0.00 450.00 \n", "2 Component_1_requirement_in_premium 0.00 -18.00 0.00 \n", "3 Component_3_requirement_in_premium 0.00 0.00 0.00 \n", "4 Component_1_requirement_in_Extra 0.00 -18.00 0.00 \n", "5 Component_2_requirement_in_Extra 0.00 0.00 -900.00 \n", "6 Component_1_availability 4500.00 20.00 0.00 \n", "7 Component_2_availability 2700.00 4.00 0.00 \n", "8 Component_3_availability 3500.00 0.00 1700.00 \n", "9 Super_minimum_production 3000.00 0.00 -0.00 \n", "10 Premium_minimum_production 3000.00 -1.20 0.00 \n", "11 Extra_minimum_production 3000.00 -6.80 0.00 \n", "\n", " Min RHS Max RHS \n", "0 -850.00 0.00 \n", "1 -450.00 Inf \n", "2 -450.00 0.00 \n", "3 -450.00 450.00 \n", "4 -450.00 0.00 \n", "5 -Inf 900.00 \n", "6 4500.00 6200.00 \n", "7 2250.00 3150.00 \n", "8 1800.00 Inf \n", "9 -Inf 3000.00 \n", "10 -0.00 3000.00 \n", "11 0.00 3000.00 " ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "total_profit = pulp.value(model.objective)\n", "display(Markdown(\"Total profit is %0.2f €\"%total_profit))\n", "\n", "display(Markdown(\"The following table shows the decision variables: \"))\n", "var_df = pd.DataFrame.from_dict(variables, orient=\"index\", \n", " columns = [\"Variables\"])\n", "var_df[\"Solution (GRB)\"] = var_df[\"Variables\"].apply(lambda item: \"{:.3f}\".format(item.solverVar.X))\n", "var_df[\"Reduced cost (GRB)\"] = var_df[\"Variables\"].apply(lambda item: \"{:.3f}\".format(item.solverVar.RC))\n", "var_df[\"Objective Coefficient (GRB)\"] = var_df[\"Variables\"].apply(lambda item: \"{:.3f}\".format(item.solverVar.Obj))\n", "var_df[\"Objective Lower bound (GRB)\"] = var_df[\"Variables\"].apply(lambda item: \"{:.3f}\".format(item.solverVar.SAObjLow) if item.solverVar.SAObjLow > -0.1 else \"-Inf\" )\n", "var_df[\"Objective Upper bound (GRB)\"] = var_df[\"Variables\"].apply(lambda item: \"{:.3f}\".format(item.solverVar.SAObjUp) if item.solverVar.SAObjUp != item.solverVar.UB else \"Inf\")\n", "\n", "\n", "display(var_df)\n", "\n", "\n", "const_dict = dict(model.constraints)\n", "con_df = pd.DataFrame.from_records(list(const_dict.items()), exclude=[\"Expression\"], columns=[\"Constraint\", \"Expression\"])\n", "con_df[\"Right Hand Side\"]=con_df[\"Constraint\"].apply(lambda item: \"{:.2f}\".format(const_dict[item].solverConstraint.RHS))\n", "con_df[\"Shadow Price\"]=con_df[\"Constraint\"].apply(lambda item: \"{:.2f}\".format(const_dict[item].solverConstraint.Pi))\n", "con_df[\"Slack\"]=con_df[\"Constraint\"].apply(lambda item: \"{:.2f}\".format(const_dict[item].solverConstraint.Slack))\n", "con_df[\"Min RHS\"]=con_df[\"Constraint\"].apply(lambda item: \"{:.2f}\".format(const_dict[item].solverConstraint.SARHSLow) if const_dict[item].solverConstraint.SARHSLow > -1e10 else \"-Inf\")\n", "con_df[\"Max RHS\"]=con_df[\"Constraint\"].apply(lambda item: \"{:.2f}\".format(const_dict[item].solverConstraint.SARHSUp) if const_dict[item].solverConstraint.SARHSUp < 1e10 else \"Inf\" )\n", "\n", "\n", "display(Markdown(\"The following table shows the constraints: \"))\n", "display(con_df)" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.7.4" }, "pycharm": { "stem_cell": { "cell_type": "raw", "source": [], "metadata": { "collapsed": false } } } }, "nbformat": 4, "nbformat_minor": 2 }